Calculation of the melting point of alkali halides by means of computer simulations 
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In this manuscript we study the liquid-solid coexistence of NaCl-type alkali halides, described by 
interaction potentials such as Tosi-Fumi (TF), Smith-Dang (SD) and Joung-Cheatham (JC), and 
compute their melting temperature (T m ) at 1 bar via three independent routes: 1) liquid/solid 
direct coexistence, 2) free-energy calculations and 3) Hamiltonian Gibbs-Duhem integration. The 
melting points obtained by the three routes are consistent with each other. The calculated T m of 
the Tosi-Fumi model of NaCl is in good agreement with the experimental value as well as with 
other numerical calculations. However, the other two models considered for NaCl, SD and JC, 
overestimate the melting temperature of NaCl by more than 200 K. We have also computed the 
melting temperature of other alkali halides using the Tosi-Fumi interaction potential and observed 
that the predictions are not always as close to the experimental values as they are for NaCl. It 
seems that there is still room for improvement in the area of force-fields for alkaline halides, given 
that so far most models are still unable to describe a simple yet important property such as the 
melting point. 



I. INTRODUCTION 

Alkali halides are inorganic compounds composed of an 
alkali metal and a halogen. The most abundant by far in 
Nature is Sodium Chloride (NaCl). NaCl in its solid form 
has a cubic structure (usually denoted as NaCl struc- 
ture) and melts at a relatively high temperature (around 
1070 K at ambient pressure) as a consequence of its high 
lattice energy. NaCl dissolves in polar solvents, such as 
water, to give ionic solutions that contain highly solvated 
anions and cations, relevant for the functioning of biolog- 
ical organisms. 

In the last century, there have been many thorough 
studies aimed at quantitatively describing the physical 
properties of alkali halides. To a first approximation al- 
kali halides are two-component mixtures of atomic anions 
and cations which interact through a spherically sym- 
metrical potential. The two unusual features of ionic 
crystals are the slow decay of the interaction potential, 
which causes long-range structural correlations, and the 
strength of the attractive cation-anion interactions. The 
simplest model for a molten salt is the restricted primi- 
tive model (RPM) , in which ions are modelled by hard 
spheres with positive and negative unit charges. The 
RPM has been studied by a number of groups 1- 7] . Back 
in 1919, Born developed a numerical model to estimate 
the energy of an ionic crystal Q. Pauling, based on 
Born's primitive model, studied the effect of the ions' 
sizes on ionic salts [t| [l(| and Mayer evaluated the role of 



polarizability and dispersive forces on alkali halides [11 1 
and proposed, together with Huggins, a generalization 
of Born's repulsive energy [l2j]. More recently, several 
ion-ion interaction potentials have been developed in or- 
der to numerically simulate alkali halides. In 1962 Tosi 
and Fumi, fitting the Huggins-Mayer dispersive energy 
to crystallographic data, proposed an empirical poten- 
tial parameterizing the repulsive part of the NaCl al- 
kali halide interactions [13] and including the dipole- 



quadrupole van der Waals term, that imitated the polar- 
ization distortion of the electronic cloud. The advantages 
and disadvantages of the Tosi-Fumi (TF) potential have 
been later underlined by Lewis and coworkers [14]. The 
model in general predicts good densities and lattice ener- 
gies of alkaline halides. However it is unable to reproduce 
dynamical features of ionic crystals such as the phonon 
dispersion curves. Also the model predicts incorrectly the 
similarity of like pair distribution functions g and g ++ 
which is absent in the experiments (TBj . The TF model 
was used for the first time in a numerical simulation by 
Adams and McDonald [la ], who obtained a remarkably 
good agreement between numerical and experimental re- 
sults. In recent years, the Tosi-Fumi potential has also 
been used to study liquid/solid phase transitions in al- 
kali halides: Valeriani et al. studied homogeneous crystal 
nucleation in molten NaCl and Chen and Zhu car- 
ried out homogeneous nucleation studies of other alkali 
halides, such as KBr [l8[ and NaBr Zykova-Timan 
et al. studied packing issues related to the interfacial 
free-energy between the liquid and solid phases (2TJ, [2l[ . 
In all these works knowledge of the melting temperature 
T m was needed. In fact, the melting point of NaCl has 
been estimated for the TF model by several groups. In 
2003, Anwar et al. determined by free-energy calcula- 
tions the melting temperature of NaCl for the Tosi-Fumi 
mode l 1221 . Later on similar results were obtained by Eike 
et al. l23l| and Mastny et al. [24| and by Zykova-Timan et 
al. [251 ] . Molecular simulations have also been performed 
to analyze experimental results (26l . [27| . In some works 
the melting temperature was calculated by direct heat- 
ing of the solid, hence the thermal instability limit rather 
that the melting temperature was calculated [Hj]. More 
recently in order to compute T m Belonoshko et al. per- 
formed two-phase simulations of NaCl and LiF [28l . [29j j . 

Besides TF, another popular model used to describe 
NaCl was proposed more recently by Smith and Dang 
(SD) [lO, [3l[ who presented an interaction potential 
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where the ion-ion interactions were Lennard- Jones like. 
This model potential has become quite popular when 
studying NaCl in water solutions, even though the prop- 
erties of its solid phase are still unknown. Similar to the 
SD potential, Joung-Cheatham (JC) presented another 
interaction potential for NaCl where the ion-ion interac- 
tions were Lennard- Jones like and proposed several NaCl 
force-fields tailored to be used in a water solution [32j . 
Somewhat surprisingly the melting point of the SD and 
JC/NaCl potentials are still unknown. 

Good model potentials for alkali halides are also use- 
ful to study ionic solutions. Simulations of alkali halides 
dissolved in water have proved useful to study thermody- 
namic mixing properties, such as the cryoscopic descent 
of the melting temperature [33j . The properties of alkali 
halides solutions have been studied at low temperatures 
in order to localize the hypothesised second critical point 
of water [34l - f3(| . The solubility of NaCl in water has also 
received certain interest [37H43I ] . However we should un- 
derline that in order to determine the solubility of a salt, 
the chemical potential of the solid should be known. 

Thus, it is relevant to quantitatively compare several 
interaction potentials in order to estimate their efficiency 
in mimicking the properties of NaCl-type alkali halides. 
To this aim, in this manuscript we evaluate the liquid- 
solid equilibrium of such systems. We have computed the 
melting temperature (at normal pressure) for three alkali 
halides, emphasizing NaCl, for three interaction models: 
the Born-Mayer-Huggins-Tosi-Fumi (TF) [IM1 13 , the 
Smith-Dang (SD) |30f and the Joung-Cheatham (JC) 
potential [32[. The three potentials are two-body and 
non-polarizable model potentials, characterized by a re- 
pulsive term, a short-range attraction and a long-range 
Coulombic interaction term. To calculate the melting 
temperature (T TO ) we used three independent routes: 1) 
liquid/solid direct coexistence 2) free-energy calculations 
and 3) Hamiltonian Gibbs-Duhem integration. We have 
found that the value obtained for the T m of the TF/NaCl 
is in good agreement with the experiment, in contrast 
with the results obtained for the SD/NaCl and JC/NaCl 
models. Therefore, we conclude that the TF/NaCl model 
is the most suitable for studies of pure NaCl. Using 
Hamiltonian Gibbs-Duhem integration we have also eval- 
uated T m for other TF/alkali halides and concluded that 
the quality of the results obtained with the TF potential 
depends on the alkali halide chosen. The TF model pro- 
vides good predictions for alkali halides that involve K + , 
Cl _ , Na + , Br - and Li + ions, whereas it performs much 
worse for alkali halides involving Rb + or F~ ions. 

The manuscript is organized as follows: we first in- 
troduce the interaction potentials of the alkali halides 
under study, i.e. Born-Mayer-Huggins-Tosi-Fumi poten- 
tial (TF), the Smith-Dang potential (SD) and the Joung- 
Cheatham (JC) potential. Then we describe the three 
simulation routes followed to compute their melting tem- 
perature: 1) the liquid/solid direct coexistence 2) the 
Einstein crystal/molecule for the solid and the thermo- 
dynamic integration for the liquid and 3) the Hamiltonian 



Gibbs-Duhem integration. To conclude we will present 
our results obtained for different alkali halides. 



II. SIMULATION METHODS 

The interaction potentials we used are two-body and 
non-polarizable model potentials, each of them charac- 
terized by a repulsive term, a short-range attractive and 
a long-range Coulomb interaction term. The TF model 
potential has the following form: 
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where r%j is the distance between two ions with charge 
<7ij, the first term is the Born-Mayer repulsive term, 

— w- and — ~w~ are the Van der Waals attractive in- 
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teraction terms and the last term corresponds to the 
Coulomb interaction. The parameters Ay, pij, Cij and 
Dij for the TF/alkali halides are given as Supplementary 
Material 

Both the SD and the JC model potentials can be writ- 
ten using the following expression: 



C/(r y ) = 4e 
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where ry is the distance between two ions with charge 
qij. The first term is Lennard- Jones-like, and its param- 
eters are given in the Supplementary Information file[45|. 
The last term in Eq. [5] corresponds to the Coulomb in- 
teraction. For the JC potential, we are going to use the 
parameters introduced to simulate NaCl in SPC/E wa- 
ter m. 

For the SD and JC the crossed interaction parame- 
ters are obtained using the Lorentz-Berthelot combining 
rules [ill |47} . It is interesting to point out that the TF 
model was developed to study ionic crystals and pure 
alkali halides in the solid phase whereas the SD was ob- 
tained to model Na + and Cl~ in water. The JC was 
fitted to model NaCl both in the solid phase and in aque- 
ous solutions. In what follows, we shall refer to ions as 
particles. 

When computing T m for the TF/NaCl interaction 
potential, we follow three independent routes: 1) liq- 
uid/solid direct coexistence, 2) free-energy calculations 
and 3)Hamiltonian Gibbs-Duhem integration. For the 
SD/NaCl and JC/NaCl and for other TF/alkali halides 
we use the first and the third route. 



A. Route 1. Liquid/solid direct coexistence 

The first route we follow to compute the melting tem- 
perature is by means of the liquid/solid direct coexis- 
tence, originally proposed in Ref. [48l452|. To start with, 
we generate an equilibrated configuration of the NaCl 
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solid phase in contact with its liquid. After having pre- 
pared the liquid-solid configuration, we run several NpT 
simulations (with anisotropic scaling, so that each side 
of the simulation box changes indepcndtly) at different 
temperatures and always at a pressure of 1 bar. Di- 
rect coexistence simulations should be performed in the 
Np z T ensemble, where z is the axis perpendicular to the 
fluid-solid interface and x and y have been chosen care- 
fully to avoid the presence of stress. However, we have 
recently shown that NpT simulations (with anisotropic 
scaling) provide proper results for sufficiently large sys- 
tems [53], [54| . Depending on the temperature, the system 
evolves towards complete freezing or melting of the sam- 
ple. Since we do not know the location of the melting 
temperature, we simulate the system in a wide range of 
temperatures and identify the melting temperature T m 
as the average of the highest temperature at which the 
liquid freezes and the lowest temperature at which the 
solid melts. A more elaborate procedure would require 
to estimate the slope of the growth rate of the solid-fluid 
interface ( which is positive at temperatures below the 
melting point and negative at temperatures above the 
melting point) and to locate the temperature at which 
the growth rate of interface is zero (i.e the melting point). 
Since the growth of the solid is a stochastic process the 
determination of the growth rate requires to accumulate 
statistics over many different trajectories making the cal- 
culations quite expensive [HI, (56[. The procedure used 
here is simpler and has prov ided reliable results for other 
systems as hard spheres [53j,|54[ or water [54j although ad- 
mittedly it yields a somewhat larger error in the estimate 
of the melting point temperature obtained from direct 
coexistence simulations (of about 5K). 



B. Route 2. Free-energy calculations 

In this route, we first compute the Helmholtz free- 
energy of both the solid and the liquid phase; next, we 
estimate the Gibbs free-energy (G) by simply adding pV . 
To compute the melting temperature we perform thermo- 
dynamic integration as a function of temperature at con- 
stant pressure to evaluate where the chemical potentials 
(p, = G/N) of both phases coincides. To estimate the 
free-energy of the bulk solid phase we use the Einstein 
crystal [53 and the Einstein molecule methods [H, H3 ■ 
Both methods are based on the calculation of the free- 
energy difference between the target solid and a refer- 
ence system at the solid equilibrium-density for the given 
thermodynamic conditions (obtained with an NpT sim- 
ulation). The reference system of the Einstein crystal 
method consists of an ideal solid whose free-energy can 
be analytically computed (an "Einstein crystal" with the 
center of mass fixed, where the inter-molecular interac- 
tions are neglected and particles are bound to their lat- 
tice equilibrium positions by a harmonic potential with 
strength Ag). The Einstein molecule method differs from 
the previous one due to the fact that we now fix only the 



position of one particle instead of the center of mass. 

Thermodynamic integration is performed in two 
steps [60j : I) we evaluate the free-energy difference (A 
Ai) between the ideal Einstein crystal and the Einstein 
crystal in which particles interact through the Hamilto- 
nian of the original solid ("interacting" Einstein crystal). 
2) Next we calculate the free-energy difference (A A2) 
between the interacting Einstein crystal and the original 
solid, by means of thermodynamic integration: U(X) — 
AJ7 soi + (l- X)(U E in-id+U S ol)), where U Em -id represents 
the energy of the interacting Einstein crystal, U so i the 
one of the original solid, and A is the coupling parameter 
that allows us to integrate from the interacting Einstein 
crystal (A = 0) to the desired solid (A = 1). The fi- 
nal expression of the Helmholtz free-energy A^ cl (T, V) 
coming from the Einstein crystal/molecule calculations 
is [60|: 

Kf l V) = A (T, V) + AAi (T, V) + AA 2 (T, V) (3) 

where Ao is the free energy of the reference system, whose 
analytical expression is slightly different in the Einstein 
crystal and Einstein molecule (see Ref. |6Cj). AAi and 
AA2 are computed in the same way in the Einstein crys- 
tal and Einstein molecule methods (the only difference 
being the choice of the point that remains fixed in the 
simulations, whether the system's center of mass or a 
reference particle's center of mass). It has been shown 
that, since the free energy of a solid is uniquely defined, 
its value does not depend on the method used to com- 
pute it and the two methodologies give exactly the same 
results [60] . It is convenient to set the thermal De Broglie 
wave length of all species to 1 A, and the internal par- 
tition functions of all species to one. These arbitrary 
choices affect the value of the free energy but does not af- 
fect phase equilibria (provided the same choice is adopted 
in all phases). 

To estimate the free-energy of the bulk liquid we use 
Hamiltonian thermodynamic integration as in Ref. [22j . 
calculating the free-energy difference between the liquid 
alkali halide and a reference Lennard- Jones (LJ) liquid 
for which the free-energy is known. Starting from an 
equilibrated NaCl liquid, we perturb the Hamiltonian of 
the system so that each ion is gradually transformed into 
a LJ atom. The path connecting both states is given by: 



U{\) = XU LJ + (I - X)U 



NaCl 



(4) 



where U LJ /JJ NaCl are the total energies of the Lennard 
Jones and NaCl fluids, respectively, and A is the coupling 
parameter (A = corresponds to NaCl whereas A=l to a 
LJ fluid). The Helmholtz free-energy of a NaCl (A^ CI ) 
is computed as: 

A%(T, V) = A% q aCl (T, V) + [~\u LJ - U NaCl ) Nm dX 



J\=o 

A™(T,F) + AA"(T,F) 



(5) 



Given that the Lennard- Jones free-energy (A^) is al- 



ready known 6I( and the integral in the Eq. [5] (AAy^) 
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can be numerically evaluated, we can determine the free- 
energy A^q Cl of the liquid alkali halide. In order to es- 
timate the integral in Eq. [5j we choose 20 values of A 
between and 1, 10 of them equally spaced from 0.000 
to 0.95 and the remaining 10 from 0.95 to 1.000, and inte- 
grate each region using the Simpson integration method. 
This choice originated from then fact that when A has a 
value close to 1.0, the integrand changes abruptly, drop- 
ping to the dispersive energy of a LJ. The Lennard- 
Jones free-energy consists of two terms: A^(T, V) — 

T, 

excess and A^' lrf (T, V) the ideal part. A^' re5 for a 
Lennard-Jones fluid has been already computed for a 
broad range of temperatures and densities in Ref. |6lLl62^ . 
The free energy of the ideal gas of a mixture of N Ma and 
Nc; is given by 



A LJM [Ti v) + A ^/ e \T, V), where A^ es (T, V) is the 
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(T, V) = N Na ln(p Na A 3 Na ) - N Na 

+ NcM Pc iA 3 cl ) - N C i 
= N[ln(p/2) - 1] (6) 



where N = Njy a + Ncz is the total number of particles 
in the system with density p = (and PNa — = 
Pci = -^t 1 = p/2) and Aj is the De Broglie thermal length 
(A; = h / ^/2TTm,ikBT), that we arbitrarily set to 1 A, to 
be consistent with our choice for the solid phase. 



C. Route 3. Hamiltonian Gibbs-Duhem integration 

The third route we follow to compute the melting tem- 
perature is by means of Hamiltonian Gibbs-Duhem ther- 
modynamic integration as in Ref. [63l - l65l ] . Starting from 
the liquid-solid coexistence point of a reference system 
(A), the Hamiltonian Gibbs-Duhem thermodynamic in- 
tegration allows one to compute the coexistence point of 
the system of interest (whose Hamiltonian is B) by re- 
solving a Clapeyron-like differential equation. In more 
detail, the Hamiltonian of the initial system (with en- 
ergy Ua), whose two-phases coexistence point is known, 
is connected to the one of the final system of interest 
(with energy f/g), via the following expression: 



U{X) = XU B + (1 - X)U A 



(7) 



where A is the coupling parameter. When two phases 
coexist (labeled as I and II): pi(T,p,X) = pn(T 7 p, A), 
being p — G/N the Gibbs free-energy per particle (the 
chemical potential of each phase). Therefore, differenti- 
ating p in both phases, we can write generalized Clapey- 
ron equations for two coexisting phases as 

v I (T,p)dp-s I (T,p)dT+ ( ^ /( J; P ' A) ) d\ = 

= v II (T,p)dp-s II (T,p)dT+ [ 8m ^' A) ] dX(8) 



where v and s are the volume and entropy per particle. 
If A is constant we recover the well known Clapeyron 
equation. If the pressure is constant we obtain the slope 
of the coexistence line in the A-T plane: 



dT = T[(dp n /d\) - (dpj/dX)} = 
dX hjj — hi 

_ T[{dp n /dX) - (dpi/dX)] 



Ah 



knowing that, at coexistence 



sii - si = 



hj j - hi 
T 



(9) 



(10) 



When a liquid coexists with a solid, (s/j — Sj) is the 
melting entropy difference As m , that can be easily com- 
puted as . Ah is obtained from the NpT simulations 

at (p,A,T) constants, whereas |j = { 9w q^ )n,p,t,\, com- 
puted with an NpT simulation at different values of A in 
each phase. 

Therefore, using Eq. [71 the generalized Clapeyron 
equation can be written as: 



dT 
~dX 



Ah 



(11) 



where ub/ua is the internal energy per particle when 
the interaction between particles is described by Ub/Ua- 
The numerical integration of the generalized Clapeyron 
equation in Eq. Ql] yields the change of the coexistence 
temperature (at constant pressure) due to the change in 
the Hamiltonian of the system, starting from the initial 
coexistence point (where interactions are described by 
Ua) to the final coexistence point (where interactions are 
described by Ub)- 



III. SIMULATION DETAILS 

When using the liquid/solid direct coexistence route 
to compute the melting temperature we choose sys- 
tems containing either 1024 [512 solid/512 liquid], 2000 
[lOOOsolid/lOOOliquid] or 4000 ions [2000 solid/2000 liq- 
uid]. When using the free-energy calculations route to 
compute the melting temperature, we simulate systems 
of 1000 ions (the unit cell of NaCl contains 8 particles). 
For the Hamiltonian Gibbs-Duhem integration route we 
simulate systems of 1000 ions. 

We simulated NaCl using the TF, the SD and the JC 
interaction potentials. For the other alkali halides con- 
sidered in this work we used only the TF potential. In 
this work we truncated the non-Coulomb part of the po- 
tential at r c =14 A and added tail corrections. We used 
Ewald sums to deal with Coulombic interactions, trun- 
cating the real part of the Ewald sums at the same cut- 
off as the non-Coulombic interactions and chose the pa- 
rameters of the Fourier part of the Ewald sums so that 
a ■ r c = 2.98 (H,^- To calculate T m for the SD and JC 
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models via direct coexistence we performed NpT molecu- 
lar dynamic simulations (MD) of systems containing 2000 
particles with the Gromacs package |68[, where we kept 
the tem per ature constant with a Nose-Hoover thermo- 
stat [o9l |7C| with a relaxation time of 2 ps, and the pres- 
sure constant at 1 bar with a Paxrinello-Rahman baro- 
stat [7l[ with a relaxation time of 2 ps. In our MD simu- 
lations, we allowed the different box lengths to fluctuate 
independently. For the TF alkali halides direct coexis- 
tence simulations we used NpT Monte Carlo simulations. 

In the free-energy calculations, we used the Einstein 
Crystal and the Einstein Molecule methods to compute 
the free-energy of the solid phase. We performed an ini- 
tial equilibration run of the solid in the NpT ensemble of 
about 10 5 Monte Carlo (MC) cycles to obtain the equi- 
librium density at the given thermodynamic conditions. 
We define a MC cycle as a translational trial-move per 
particle and a trial-volume change. For the thermody- 
namic integration in the NVT ensemble we carried out 
2xl0 4 equilibration and 8xl0 4 production cycles for every 
value of A and simulated 20 values of A per thermody- 
namic state. We also used thermodynamic integration to 
compute the free-energy of the liquid phase. We carried 
out 8xl0 4 equilibration and 18xl0 4 production cycles for 
every value of A and simulated 21 values of A per thermo- 
dynamic state. Free-energy calculations were performed 
at 1083 K and 1 bar. To obtain the equilibrium den- 
sities at these thermodynamic conditions, we run NpT 
MC simulations of the liquid and solid phases. Once 
equilibrated, we used those densities in the free energy 
calculations. 

When performing the thermodynamic integration to 
compute the free-energy of the liquid phase (route 2), 
we tested the dependence of our results on the choice of 
the reference system by performing the thermodynamic 
integration to two Lennard- Jones models with different 
parameters: the parameters set LJ1 used by Anwar et 
al. [22j and another set of parameters denoted as LJ2. 
Both sets of parameters are represented in Table |U 



e/k B [LJ1] 


a 


P* 




e/k s [ LJ2] 


a 


P* 




537.01 


2.32 


0.3766 


2.02 


358.00 


3.00 


0.8143 


3.03 



TABLE I. LJ1 and LJ2 sets of parameter with units 
e/k B ([A"]) and a{[A}). p*=pa 3 and T*=fc s T/e. 

The values of p* and T* shown in Table U are ob- 
tained by scaling the density and temperature of the liq- 
uid phase of the TF/NaCl at 1083 K and 1 bar to LJ 
reduced units. The free-energy should be independent of 
the choice of the LJ reference system. Notice that the 
free-energy of the LJ system given by the Nezbeda equa- 
tion of state (EOS) |61| has a lower error when using LJ2 
(at p*=0.8143 and T*=3.03) rather than at p*=0.3766 
and T*=2.02 (when using LJ1) due to the proximity of 
the LJ fluid critical point. 

When using the third route to compute the melting 



temperature, we integrated the Hamiltonian in Eq. [7] 
from the TF/NaCl to the SD and JC potentials and from 
the TF/NaCl potential to other alkali halides potentials 
parameterized using TF. In all cases, we simulated 5 val- 
ues of A per Hamiltonian Gibbs-Duhem integration. 



IV. RESULTS 

Let us start by presenting the results for the melting 
point of NaCl for the models and routes considered in 
this work. 



A. Route 1. Liquid/solid direct coexistence 

Using the liquid/solid direct coexistence technique fill 
|49| we determined the melting temperature at 1 bar not 
only of the TF/NaCl but also of the other NaCl poten- 
tials (the Smith-Dang and Joung-Cheatham). After hav- 
ing prepared the equilibrated liquid-solid configuration 
we performed several anisotropic NpT Monte Carlo sim- 
ulations (in the case of TF/NaCl) for two different sys- 
tem sizes, N=1024 and N=4000, to analyze the finite-size 
effects. For SD/NaCl and JC/NaCl we used molecular 
dynamics simulations with an anisotropic barostat in sys- 
tems containing 2000 particles. In Fig. Q] we plot the 
time evolution of the total energy of the TF/NaCl sys- 
tem equilibrated at 1 bar and at different temperatures 
for the two system sizes. 
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FIG. 1. Total energy (in kcal per ion mol) versus MC cycles 
for the Tosi-Fumi NaCl at T= 1085 K, 1084 K, 1083 K and 
1082 K for the 1024 particles system (left-hand side) and T= 
1100 K, 1082 K and 1075 K for the 4000 particles system. 

After an equilibration interval of about 10 4 MC cy- 
cles (where the energy stays constant), we observed that 
when the temperature is below melting the energy de- 
creases until it reaches a plateau with a sudden change of 
slope, corresponding to the situation where the liquid has 
fully crystallized. When the temperature is above melt- 
ing, the energy increases until it reaches a plateau and 
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stays constant: at this stage, the solid has completely 
melted. The results obtained for the 1024 particles sys- 
tem (left-hand side of Fig. [J) show that when T=1083 K 
the liquid crystallizes, whereas when T > 1084 K the solid 
melts: therefore, the estimated melting temperature for 
the 1024 particles system is T m =1083(5) K. The results 
obtained for the 4000 particles system (right-hand side 
of Fig. [IJ show that when T < 1075 K the liquid crys- 
tallizes, whereas when T > 1100 K the solid melts. At 
T=1082K the energy remains approximately constant. 
From this we conclude that for the 4000 particles system 
T TO =1082(13) K. Thus finite size effects on the melting 
point of NaCl seems to be small once the system has at 
least 1000 particles. 

Hence, when computing T m for the remaining interac- 
tion potentials, we chose a large enough system with 2000 
particles. In Fig. [5J we present the time evolution of the 
total energy of the NaCl/SD and NaCl/JC 2000 parti- 
cles systems, at 1 bar and at different temperatures. The 




t [ps] t [ps] 

FIG. 2. Total energy (in kcal per ion mol) versus time 
(in picoseconds) for the NaCl/SD system (left-hand side) at 
T=1329 K, 1328 K and 1327 K, and for the SD/NaCl system 
(right-hand side) at T= 1284 K, 1285 K and 1286 K. Note 
the different y-axis of both plots. 

results obtained for the 2000 particles system show that 
the estimated melting temperature is T m = 1327(5) K for 
the NaCl/SD and T TO =1285(5) K for the NaCl/JC, re- 
spectively 

From these results, we can already conclude that the 
potential that gives the melting temperature closest to 
the experimental one (1074 K) is the Tosi-Fumi. For this 
reason the TF/NaCl model is the most suitable for simu- 
lations of pure NaCl. This conclusion is further confirmed 
when calculating the melting curve for both TF/NaCl 
and JC/NaCl potentials. Using Gibbs Duhem integra- 
tion [72J we calculated the p — T melting curve of the 
TF/NaCl and the JC/NaCl, presented in Fig.H 

Concerning the results of the TF/NaCl, we observe 
that at low pressure our calculations recover the experi- 
mental slope of the melting curve = 30.6(5) bar/K), 




2000 



T(K) 

FIG. 3. Melting curve for the TF/NaCl (open diamonds) 
and JC/NaCl (open circles). Experimental results are from 
Ref. [2q ] (open squares) and Ref. [7l| (open triangles). 

in good agreement with previous calculations [22|, [24j], 
whereas at higher pressures the slope of the melting line 
is lower than the experimental one. On the other hand, 
concerning the results for the JC/NaCl, we observe that 
not only the melting temperature, but also the slope of 
the melting curve does not reproduce the experimental 
one. The same conclusions can be drawn for the melting 
enthalpy (Ah m ). The calculated melting enthalpies for 
the TF/NaCl, SD/NaCl and JC/NaCl are 3.36 kcal/mol, 
4.8 kcal/mol and 4.6 kcal/mol, respectively, that com- 
pared to the experimental value (3.35 kcal/mol), confirms 
the better performance of the TF/NaCl model with re- 
spect to the other models. 

B. Route 2. Free-energy calculations 

According to the second route, we first computed the 
Helmholtz free-energy of the solid and liquid phase of the 
TF/NaCl, and then estimated the Gibbs free-energy of 
each phase (G) by adding the pV term. After that, we 
performed thermodynamic integration of G in the (p,T) 
plane at constant pressure as in Eq. 1121 

G(T 2 , P ) = G(Ti,p) _ [ T * H(T) 

Nk B T 2 Nk B T x J Tl Nk B T 2 [ ' 

where the enthalpy H(T) = U(T)+pV(T) can be easily 
obtained in NpT simulations at each temperature. The 
melting point is defined as the state point where the two 
phases have the same chemical potentials (fi = G/N). 

We computed the free-energy of the NaCl solid phase 
at 1083 K and 1 bar with Einstein crystal (EC) [53| and 
Einstein molecule (EM) [58| algorithms. Our results are 
summarized in Table UH where we present free-energies 
in Nk B T units. To determine the normal melting point 
any temperature could have been selected, the choice of 
1083 K is convenient from a practical point of view since 
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it is expected to be close to the T m of the model, so that 
to reduce the contribution of the second term on the right 
side of the Eq. QH 





r[K] 


N 


P (N/A 3 ) 




Aa[NksT] 


A Ai [NksT] 


A A 2 [Nk fl T] 


Af„f [NksT] 


TF/EC 


1083 


1000 


0.03856 


500 


7.583 


-42.73 


-6.33 


-41.181(9) 


TP/EM 


1083 


1000 


0.03856 


500 


7.594 


-42.73 


-6.34 


-41.477(9) 



TABLE II. Free energy from the Einstein Crystal/Einstein 
Molecule for the TF/NaCl solid phase at 1 bar and different 
temperatures (T) and system size (N being the total number 
of particles), p the number density (N/V), Ae the spring 
constant, Ao, AAi, AA2 are the terms in Eq. 0and the free 
energy of the solid A^°[ CI is represented in the last column 
and corresponds to Eq[3] 

In Table UH we observe a perfect agreement between 
the free-energy computed with the Einstein Crystal and 
the one compute with the Einstein Molecule at the same 
pressure, temperature and system size (1000 particles). 

Having computed the chemical potential of the solid 
phase, we calculated the Helmholtz free-energy of the 
liquid phase using two Lennard- Jones reference systems 
LJ1 and LJ2 in order to study the uncertainties associ- 
ated to the choice of the reference system. Results are 
summarized in Table Hill 





AA^[Nk B T] 


A LJ M[mBT] 


^Ha' re8 [NksT] 


A~° c, [Nk B T] 


TF/LJl 


35.95 


-5.19 


-0.327 


-41.47(2) 


TF/LJ2 


37.12 


-5.19 


0.837 


-41.48(2) 




FIG. 4. (U LJ - U NaCl ) N ,v,T,x from Eq. [5] at 1083 K and 1 
bar. The two curves correspond to a different set of Lennard- 
Jones parameters for the reference system: LJ1 (open circles) 
and LJ2 (open squares). 



MC simulations to compute enthalpy H(T) and to esti- 
mate the equation of state of the liquid and solid phases 
(where a typical MC run consists of 3xl0 4 equilibration 
and 7x10 production cycles). The chemical potential 
of one phase is fi = G/N and the melting temperature 
is given by the point at which the two phases have the 
same chemical potentials. Our results for the TF/NaCl 
are presented in Table HVl 



TABLE III. TF/NaCl free energy of the liquid phase (A^ ct ) 
as in Eq. [5] at 1 bar and T=1083 K. The results presented 
in the table refer to a system with N = 1000 particles and 
number density of p — 0.03016 A~ 3 . AA LJ is the integral in 

A L.l,id 

Eq-U jjjfpr = ln(p/2) - 1.00. Af/' res is obtained from the 
EOS of Ref. [H. 

As shown in Table [TTT1 the value of the integral ( AA^ ) 
depends on the choice of the LJ parameters of the ref- 
erence LJ system. In Fig. @] the integrand of Eq. [5] is 
shown when carrying out the thermodynamic integration 
from the liquid TF/NaCl to both LJ reference systems 
LJ1 and LJ2. It is relevant to stress that, due to the 
abrupt change of (U LJ — JJ NaCl ) n.v.t.x for values of A 
close to 1, many points were used in the integration be- 
tween A = 0.95 and A = 1.00. 

To calculate the residual free energy of the reference 
LJ fluid at our thermodynamic conditions, we used the 
equation of state (EOS) for the LJ system proposed by 
Nezbeda et al. [6l|, 111]- Independently on the chosen 
reference system (whether LJ1 or LJ2), the two free- 
energies A^ aC7 [Nk B T] coincide (last column in Table 

En}. 

After having computed the Helmholtz free-energy, we 
estimated the Gibbs free-energy (G) by adding the pV 
term and performed thermodynamic integration of G 
at constant pressure (see Eq. [P2)) where we used NpT 





TmpC] 




T m [K] 


TF/LJl 


1083(3) 


TF/LJ2 


1084(3) 



TABLE IV. Melting temperature at 1 bar of the TF/NaCl 
model. The results presented in the table refer to a system 
with N — 1000 particles. The parameters for the interaction 
potentials of the reference systems (LJ1 and L J2) are given 
in Table |U 

Both results are in perfect agreement with the T m cal- 
culated by direct coexistence. Other EOS could be used 
to calculate the residual free energy of the reference LJ 
fluid, such as the one proposed by Johnson et al. [75$ . 
When the LJ residual contribution is taken from John- 
son et al, the melting temperature turns out to be about 
6 K higher than when it is taken from Nezbeda et al. 
However, it is likely that the EOS proposed by Nezbeda 
is slightly more accurate than that by Johnson et al. [6l[ . 
In any case, the differences are small. 



C. Route 3. Hamiltonian Gibbs-Duhem integration 

The third route we followed to compute the melting 
temperature of several alkali halides at 1 bar is by means 
of the Hamiltonian Gibbs-Duhem thermodynamic inte- 
gration. We used the TF/NaCl as the reference Hamilto- 
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nian (Ua in Eq.[7]) and integrated the generalized Clapey- 
ron equation (Eq.QT]) using a 1000-particle system. From 
the A = 1 point of the Hamiltonian Gibbs-Duhem in- 
tegration, we could compute the melting entropy A5 m 
(being AS m = NAs m , see Eq. [T0|) and the melting 
enthalpy (AH m = NAh m ) at coexistence. In the Ta- 
ble |V] we present our results for the melting temper- 
ature, AS m , reporting also the experimental values of 
both for each alkali halide, and for AH m . We checked 
our Hamiltonian Gibbs Duhem calculations by comput- 
ing the melting temperature of TF/KF at 1 bar by means 
of liquid/solid direct coexistence. Figure [5] represents the 
time evolution of the total energy of the TF/KF 2000- 
particle systems, equilibrated at 1 bar and at different 
temperatures. From these results we estimate T m to be 
860(5) K, that corroborates the result obtained with the 
Hamiltonian Gibbs Duhem integration in Table |V] (i.e. 
859(15)K). 



Alkali halide 


T 


m 


AS m 




AH m 




NaCl 


1082(13) 


1074 


3.10 


3.12 


3.36 


3.35 


KF 


859(15) 


1131 


2.77 


2.98 


2.39 


3.37 


KBr 


1043(15) 


1003 


3.14 


3.03 


3.29 


3.04 


KC1 


1039(15) 


1049 


3.12 


3.04 


3.25 


3.19 


RbBr 


1047(15) 


955 


3.34 


2.88 


3.51 


2.75 


RbCl 


1092(15) 


988 


3.18 


2.88 


3.48 


2.85 


RbF 


996(15) 


1048 


3.03 


2.88 


3.03 


3.02 


NaF 


611(15) 


1266 


2.10 


3.09 


1.29 


3.91 


NaBr 


1022(15) 


1028 


3.10 


3.06 


3.18 


3.70 


LiCl 


780(15) 


887 


2.46 


2.70 


1.93 


2.39 


LiF 


1010(15) 


1118 


2.93 


2.88 


2.97 


3.22 



TABLE V. Melting temperatures, AS m and AH m at 1 bar 
for several alkali halides simulated using the TF interac- 
tion potential with units T m ([.K]), AS m ([cal / (Kmol)]) and 
AHm, ( \kcal/mol] ) . The experimental values of AS m are from 
Ref. [76H78I]. The bold fonts represent alkali halides for which 
there are good agreement with the experimental melting tem- 
peratures. 

While the Tosi-Fumi potential predicts the NaCl melt- 
ing temperature in good agreement with its experimental 
value, the results for the other alkali halides are not as 
accurate. The predictions for KC1 and NaBr are reason- 
able (differing only about 15 K and 10 K from the experi- 
ments). Whereas for alkali halides involving Rb or F ions 
(such as KF, RbBr, RbCl, RbF, LiF and NaF) the cal- 
culated values strongly differ from the experimental ones 
(with deviations of up to 100 K). From this we conclude 
that the TF model potential, being non-polarizable, fails 
when describing alkali halides that involve big cations 
(such as Rb + ) and small anions (such as F~). In general, 
when T m of a given TF/alkali halide is lower/higher than 
the experimental value, also AH m is smaller/larger than 
the experiments, so that the predicted AS m is in reason- 
able agreement with the experimental value for most of 
the studied alkali halides. 




MC cycles 

FIG. 5. Total energy versus MC cycles for TF/KF at T = 
865 K, 860 K, 855 K and 850 K (from top to bottom). 

V. DISCUSSIONS 

We now compare our data for the melting tempera- 
ture of NaCl at 1 bar with other values taken from the 
literature (see Table |VI*|) . 





T m [K] 


N 


Pliq 


Psol 


TF/ Anwar [22] 


1064(14) 


512 






TF/ Anwar* [22] 


1084(14) 


512 






TF/Zykova-Timan [25] 


1066(20) 


5760 






TF/Mastny [24] 


1050(3) 


512 






TF/Eike [23] 


1089(8) 


512 






TF/An [79] 


1063(13) 


4096 






TF/this work 


1083(5) 


1024 


1.465 


1.876 


TF/this work 


1082(13) 


4000 


1.465 


1.876 


SD/this work 


1327(10) 


2000 


1.216 


1.668 


JC/this work 


1286(10) 


2000 


1.283 


1.746 


Experiments/ Janz [801 


1074 




1.54 


1.98 



TABLE VI. Melting temperature of the NaCl computed using 
TF, SD and JC interaction potentials at 1 bar. N represents 
de number of ions of each study. pu q and p so i are the liquid 
and solid densities in gcm~ 3 , respectively. The * represents 
the recalculated T m at 1 bar using the values given by Anwar 
(p CO e:r=-300 bar at 1074 K and |£ = 30 bar K' 1 ). 

Concerning the calculations of the melting tempera- 
ture with the Tosi-Fumi, we observe that T m calculated 
in this work is in good agreement (within the error bar) 
with the one reported by Anwar et al. [HJ . They calcu- 
lated the coexistence pressure at 1074 K and —300 bar. 
Then, they evaluated the slope of the coexistence curve, 
■J* = 30 bar K -1 , and recalculated the melting temper- 
ature at 1 bar obtaining 1064(14) K (nonetheless, using 
their values, we have obtained the NaCl melting tem- 
perature at 1084(14) K, in perfect agreement with our 
results). Zykova-Timan et al. [25j j computed the melting 
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temperature via liquid/solid direct coexistence and ob- 
tained a value for T m in good agreement with the one 
obtained in this work (taking into account that in their 
NVE runs they obtained the melting temperature for a 
pressure of about +- 500 bar and since the slope of the 
melting curve is of about 30K/bar the value at p = lbar 
could be modified by up to 17K ). In general, the agree- 
ment with data in the literature is satisfactory. The TF 
slightly underestimates the experimental values of the 
coexistence densities, whereas both SD and JC underes- 
timate them considerably, in agreement with the results 
of Alejandre and Hansen for the SD model [8lj |. 

When computing the melting temperature of the 
TF/NaCl model potential via liquid/solid direct coexis- 
tence and free-energy calculations we have studied finite- 
size effects and observed that the value of the melting 
temperature did not change much for systems having 
more than 1000 particles. Thus we take T m =1082(13) K 
(i.e the melting point of the 4000 particle system ) as 
our estimate of the melting point of the TF/NaCl model 
potential in the thermodynamic limit ( where the error 
bar now includes the possible contribution to the error of 
the extrapolation to infinite size). When calculating T m 
via Hamiltonian Gibbs-Duhem integration, we have used 
the value of T m =1082(13) K as the initial point for the 
Hamiltonian Gibbs Duhem integration at coexistence. 

Concerning the calculation of the melting temperature 
with the Smith-Dang and Joung-Cheatham potentials, 
we computed these values via two independent routes 
(liquid/solid direct coexistence and Hamiltonian Gibbs- 
Duhem integration) and obtained the same results. The 
T m of both models is reported here. From our calcu- 
lations we conclude that these models overestimate the 
melting temperature of NaCl, being about 200-250 K 
higher than the experimental value of 1074 K. Therefore, 
it is clear that the melting temperature that most resem- 
bles the experimental one at 1 bar is the one calculated 
using the TF/NaCl. Although JC/NaCl or SD/NaCl 
would work for NaCl solutions in water, they seem to 
be unsuitable for simulations of pure NaCl. 



VI. CONCLUSIONS 

In this manuscript, we have computed the melting tem- 
perature at 1 bar for different NaCl-type alkali halides. 
When computing T m for the Tosi Fumi NaCl interaction 
potential, we have followed three independent routes: 1) 
liquid/solid direct coexistence, 2) free-energy calculations 
and 3) Hamiltonian Gibbs-Duhem integration. For the 
Smith Dang /NaCl and Joung Cheatham NaCl, we have 
used the first and third route, whereas for other Tosi 
Fumi alkali halides we have applied only the third route. 

The results obtained for the T m of TF/NaCl are in 
good agreement with other numerical [2(J, [22[ and with 
experimental results [8(3], giving T m = 1082(13) K at 
1 bar. When computing T m for the SD/NaCl and 
JC/NaCl, we find a perfect agreement between the cal- 



culations obtained via liquid/solid direct coexistence and 
Hamiltonian Gibbs-Duhem integration. However, both 
models overestimate the melting temperature of NaCl 
by more than 200 K. We have also determined the melt- 
ing curve for the Tosi-Fumi and Joung-Cheatham models 
and found that the Tosi-Fumi correctly predicts the be- 
havior of the curve (Jfr) at low pressures, but does not 
capture the experimental behavior when the pressure in- 
creases. Therefore, we conclude that the SD/JC models 
are unable to reproduce the properties of pure NaCl. 

We have also computed the melting temperature of 
other alkali halides using the Tosi-Fumi interaction po- 
tential and observed that this model gives good predic- 
tions for NaCl, NaBr, KC1 and KBr, whereas for the other 
alkali halides the predictions are not as good, especially 
when it concerns Rb + and F~ ions. The reason for this 
probably being that the Tosi-Fumi interaction potential 
is not polarizable and cannot capture the highly polariz- 
able character of these ions. Neglecting the polarization 
terms causes an incorrect description of these salts. In 
the past polarization effects were normally incorporated 
into the simulations of ionic systems via the shell model 
f82j ]. More elaborate models developed recently employ 
potentials which include the polarization effects using 
either multipole expansions 15] or the distortable ion 
model [83]. The TF/alkali halide potentials have a seri- 
ous transferability problem: the same ions present differ- 
ent potential parameters depending on the alkali halide 
(i.e. the Na-Na interaction is different in NaCl and NaF). 
In the case of LJ-like models, it seems that it is not pos- 
sible to describe NaCl accurately with a model consisting 
of point charges and a Lennard- Jones interaction site for 
each ion where Lorentz-Berthelot rules are used to de- 
scribe the interaction between cations and anions. The 
same conclusion was obtained by Cavallari et al. who 
showed how the ion-ion interaction reconstructed from 
standard LB rules fail to correctly account for the struc- 
ture of the concentrated solutions model with ab-initio 
calculations. An interesting possibility would be incor- 
porating deviations to the Lorentz-Berthelot combining 
rules to obtain the crossed interactions between the ions 
or adjusting these interactions as in the TF model, al- 
though probably the use of a Van der Waals r -6 term 
would not capture the underlying physics. Thus, it seems 
there is still room for improvement in the area of alkali- 
halides salt models. 
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